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ABSTRACT 


The  Schur  algorithm  and  its  times-domain  counterpart,  the  fast 
Cholseky  recursions,  are  some  efficient  signal  processing  algorithms 
which  are  well  adapted  to  the  study  of  inverse  scattering  problems. 

These  algorithms  use  a  layer  stripping  approach  to  reconstruct  a  loss¬ 
less  scattering  medium  described  by  symmetric  two-component  wave  equations 
which  model  the  interaction  of  right  and  left  propagating  waves.  In  this 
paper,  the  Schur  and  fast  Cholesky  recursions  are  presented  and  are  used 
to  study  several  inverse  problems  such  as  the  reconstruction  of  nonuniform 
lossless  transmission  lines,  the  inverse  problem  for  a  layered  acoustic 
medium,  and  the  linear  least-squares  estimation  of  stationary  stochastic 
processes.  The  inverse  scattering  problem  for  asymmetric  two-component 
wave  equations  corresponding  to  lossy  media  is  also  examined  and  solved  by 
using  two  coupled  sets  of  Schur  recursions.  This  procedure  is  then  applied 
to  the  inverse  problem  for  lossy  transmission  lines. 
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I .  Introduction 

The  Schur  algorithm  [1]  -  [2]  is  a  fast  algorithm  well-suited  to  high¬ 
speed  data  processing.  This  algorithm  is  obtained  by  using  a  layer  stripping 
procedure  to  reconstruct  a  lossless  scattering  medixim  described  by  symmetric 
two-component  wave  equations.  Two-component  wave  equations  describe  many 
real-world  physical  phenomena,  and  the  purpose  of  this  paper  is  to  examine 
several  inverse  problems  which  can  be  solved  efficiently  by  using  the  Schur 
algorithm.  The  problems  that  we  will  consider  include  the  reconstruction  of 
nonuniform  lossless  transmission  lines,  the  inverse  problem  for  an  acoustic 
layered  medium  and  the  linear  least-squares  estimation  of  stationary  stochastic 
processes. 

In  addition,  we  will  consider  the  case  when  the  scattering  medium  that 
we  want  to  reconstruct  is  lossy  and  is  described  by  asymmetric  two- component 
wave  equations.  It  will  be  shown  that  in  this  case  the  medium  can  be  recon¬ 
structed  by  using  coupled  Schur  recursions,  and  this  procedure  will  be 
illustrated  for  the  case  of  nonuniform  lossy  transmission  lines. 

The  Schur  algorithm  has  a  breadth  of  applications  which  is  nothing 
short  of  astonishing.  It  is  used  for  example  in  the  Schur-Cohn  test  [3] 
for  checking  the  stability  of  discrete-time  polynomials ,  and  in  the  Darlington 
procedure  [3] ,  [4]  for  synthesizing  an  impedance  function  by  a  cascade  of 
elementary  lossless  sections  of  degree  one  terminated  by  a  resistor.  The 


fast  Cholesky  recursions,  which  constitute  the  time-domain  version  of  the 
Schur  algorithm,  may  be  used  to  obtain  a  lower  triangular-diagonal-upper 
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triangular  (LDU)  factorization  of  a  Toeplitz  matrix  [5]  -  [7] .  The  continuous 
parameter  version  of  this  algorithm,  which  will  be  employed  throughout 
this  paper,  similarly  performs  a  causal- anticausal  factorization  of  a  Toeplitz 
operator  [8] . 

This  paper  is  organized  as  follows.  In  Section  II,  we  consider  the  problem 
of  reconstructing  a  lossless  scattering  medi\m  described  by  symmetric  two- 
component  wave  equations  and  we  show  that  the  medium  can  be  reconstructed 
layer  by  layer  by  using  the  fast  Cholesky  recursions  or  the  Schur  algorithm. 

This  procedure  is  then  used  in  Section  III  to  solve  the  inverse  problem  for 
nonuniform  lossless  transmission  lines,  and  scattering  concepts  such  as 
those  of  impedance,  reflection  coefficient,  right  and  left  propagating  waves 
are  illustrated.  In  Section  IV  the  inverse  problem  for  an  acoustic  layered 
medium  is  examined,  and  the  Schur  algorithm  (which  is  known  in  this  context 
as  the  dynamic  deconvolution  algorithm  [9])  is  used  to  reconstruct  the 
medium  parameters  from  scattering  data.  The  data  are  obtained  by  probing 
the  medium  with  plane  waves,  first  at  normal  incidence  and  then  at  non¬ 
normal  incidence.  In  Section  V,  the  problem  of  linear  least-squares 
estimation  of  a  stationary  stochastic  process  in  white  noise  is  discussed.  It  is 
shown  that  the  forwards  and  backwards  estimation  residuals  satisfy  two  component 

wave  equations  and  that  the  lattice  estimation  filter  obtained  by  discretizing 
these  equations  can  be  constructed  from  the  spectral  function  of  the  observed 
process  by  using  the  Schur  algorithm.  In  Section  VI,  the  inverse  problem 
for  a  lossy  scattering  medium  is  considered,  and  is  solved  by  using  two 
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coupled  Schur  recursions.  As  an  example  of  this  result,  the  inverse  problem 
for  lossy  nonuniform  transmission  lines  is  examined.  Section  VII  contains 
some  conclusions  and  suggests  topics  for  future  research. 

Notation 

Unless  specifically  identified  otherwise,  all  variables  are  scalar. 

Fourier  transforms  will  be  designated  by  a  carat,  e.g.  x,  and  partial 

d^f 

derivatives  by  subscripts;  =  f  .  Dependent  variables  will  generally 

CliXClTl 

be  omitted  for  brevity. 
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II .  The  Continuous-Time  Schvir  and  Fast  Cholesky  Algorithms 

In  this  section,  we  quickly  review  the  continuous- time  Schur  and  fast 
Cholesky  algorithms.  No  derivations  will  be  attempted  since  all  results 
can  be  found  in  the  given  references.  We  consider  a  lossless  scattering 
medium  described  by  the  symmetric  two-component  wave  equations 

=  -r(x)q(x,t)  (la) 

-  <1^  =  -r(x)p(x,t)  (lb) 

which  constitute  a  special  case  of  the  equations  discussed  in  [10]  -  [11] . 
Here  r(x)  is  the  reflectivity  function  and  p(x,t)  and  q(x,t)  are  the 
rightward  and  leftward  propagating  waves  in  the  medixm  at  point  x  and 
time  t.  Note  that  if  r(x)  s  0  over  some  interval,  then 

p  =  p(x-t)  q  =  q(x+t)  (2) 


over  the  interval,  so  that  p  and  q  correspond  effectively  to  waves 
propagating  rightward  and  leftward  with  unit  velocity. 

In  the  following,  it  will  be  assumed  that  r (x)  E  0  for  x<0  and  that 
r  e  L^[0,oo)  ,  so  that  for  x<0  and  as  x-x»,  p(x,t)  and  q(x,t)  have  the  form  (2)  , 

The  Scattering  Matrix 


By  taking  the  Fourier  transform  of  (1) ,  we  obtain 


d 

dx 


-jw  -r(x) 

-r(x)  j(jj 


(3) 


and  a  simple  discretization  of  x  in  (3)  gives  the  elementary  scattering 
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section  described  in  Figure  1.  This  figure  shows  that  r(x)A  is  the  fraction 
of  the  rightgoing  wave  p  which  is  reflected  by  a  section  of  thickness 
A  at  depth  x  inside  the  mediimi.  The  discrete  ladder  structure  displayed 
by  Figxire  1  has  been  used  to  design  signal  processing  architectures  for 
^eech  processing  [12] ,  digital  wave-filter  synthesis  [13] ,  spectral 
estimation  [14]  and  linear  estimation  [4],  [15]  -  [16], 

The  elementary  scattering  layers  of  Figure  1  can  be  conposed  by  using 
the  rules  of  conposition  for  scattering  layers  described  in  Redheffer 
[17] .  The  resulting  aggregate  mediimi  is  described  by  the  scattering 
matrix 


3(0)) 


T^  (W) 

Lt 

R^(a3) 

Li 

^j^(w) 

(4) 


which  relates  the  incoming  and  outgoing  waves  appearing  in  Figures  2a  and 
2b.  In  Figure  2a,  the  itiediTom  is  probed  from  the  left  by  a  rightward 
propagating  wave  e  and  R^(u)e^^  and  'f^(aj)e  are  respectively 
the  reflected  and  transmitted  waves.  Figure  2b  corresponds  to  the  case 
when  the  mediiim  is  probed  from  the  right.  More  generally,  for  arbitrary 
waves  p(x,(jo)  and  q(x,a)) 


p(x,u)  =  p  {(jo)e  q(x,a))  =  q  (aj)e^^ 

Ij  Li 


(5a) 


for  x<0 ,  and 

p(x,ca)  =  Pj^(oo)e  q(x,a))  =  qj^((ja)e 


(5b) 


and  x-x»,  and 
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“  ^  “1 

Pjj(to) 

=  S((jO) 

PL(a)) 

q^  (U) 

Xj 

q^((.) 

(6) 


expresses  the  outgoing  waves  (p  ,  q  )  in  function  of  the  incoming  waves 

K  Xj 

,/S  » 

'Pl'  V- 

If 


(x,(iO) 


p^{x,a)) 


(x,a)) 


i  =  1,2 


(7) 


are  two  arbitrary  solutions  of  (3),  and  if  S  =  diag(l,-l),  the  system  (3) 
has  the  property  that 


(a^(x,oj)  E  a  (x,u)))  =  0 
dx  1  2 


(8) 


and 

W(a^  (XjW)  ,  a^Cx,^))  =0,  (9) 

where  H  denotes  the  Hermitian  transpose,  and  where 

W(a  ,a  )  » 

is  the  Wronskian  of  a^  and  a^-  The  relation  (8)  expresses  the  fact 
that  the  medium  is  lossless  and  it  can  be  used  to  show  that  the  matrix 
S((jj)  is  unitary,  i.e. 

S^(w)S{a))  =  I  .  (11) 


This  property  is  valid  only  when  the  two  component  system  (3)  is  symmetric, 
whereas  the  identity  (9)  holds  even  for  asymmetric  two-component  systems 
of  the  type  that  will  be  considered  in  Section  VI . 

A  consequence  of  (9)  is  that 

T^  (OJ)  =  Tj^  ((jO)  . 


(12) 
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Physically,  this  means  that  the  transmission  loss  through  the  system  is  the 
same  going  in  either  direction.  It  can  also  be  shown  (see  [10])  that  the 
assumption  that  r  6  L  [0,°°)  implies  that  T  (co)  has  no  poles  in  the  lower 
half -plane,  so  that  the  system  (3)  has  no  bound  states.  By  combining  this 
observation  with  (11)  and  (12) ,  we  can  conclude  [18] ,  [19]  that  the  entries 
of  S(a3)  can  all  be  computed  from  either  R  (w)  or  R  (w)  .  This  property  is 
very  important,  since  in  some  inverse  scattering  applications  such  as  the 
inverse  seismic  problem,  we  have  access  to  only  one  side  of  the  scattering 
medium.  Note  that  our  sign  convention  for  the  Fourier  transform  is  the 
opposite  of  that  of  [18]  -  [19] ,  which  explains  why  we  use  the  lower  half-plane 
to  study  the  properties  of  S(a)),  instead  of  the  upper  half-plane  in  [18]  -  [19]. 

Fast  Cholesky  Recursions 

To  obtain  the  fast  Cholesky  recursions,  we  assume  that  the  medi\mi  is 
quiescent  at  t=0,  and  that  is  is  probed  from  the  left  by  a  known  rightward 
propagating  wave 

p(.0,t)  =  6  (t)  +  p(0,t)u(t)  (13) 

which  is  incident  on  the  mediiim  at  t=0.  Here  5(-)  denotes  the  Dirac  delta 
function  and 

1  for  t^O 

uCt)  =  (14) 

0  for  t<0 

is  the  unit  step  function.  Note  that  the  main  feature  of  p(0,t)  is  that 
it  contains  a  leading  impulse  which  is  used  as  a  tag  indicating  the 
wavefront  of  the  probing  wave.  The  measured  data  is  the  reflected  wave 
q(0,t)  =  q(0,t)u(t) 

recorded  at  x=0.  In  the  special  case  when  p(0,t)  H  0,  q(0,t)  =  R  (t) 

Xj 


(15) 
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is  the  inpulse  response  of  the  scattering  medium  and  its  Povirier  transform 
is  the  left  reflection  coefficient.  Note  that  ^  can  also  be 
measxired  by  sending  into  the  medim  sinusoidal  waveforms  at  various 
frequencies  and  measuring  the  magnitude  and  phase  shift  of  the  reflected 
sinusoidal  wave.  In  the  following,  for  convenience  we  will  omit  the 
subscript  L  of  R  (t)  and  R  (to)  . 

Xj  Xj 

Since  the  medium  is  causal  and  originally  at  rest,  the  waves  p(x,t) 
and  q(x,t)  inside  the  medium  must  have  the  form 


p(x,t)  =  6(t-x)  +  p(x,t)u{t-x) 


(16) 


q(x,t)  =  q(x,t)u(t-x) 

where  p(x,t)  and  q(x,t)  are  smooth  functions.  By  substituting  (16)  inside 
(1) ,  and  identifying  coefficients  of  the  iipulse  6(t-x)  on  both  sides 
of  (lb) ,  we  find  that 


r{x)  =  2q(x,x)  (17) 

and 

=  -r(x)q(x,t)  (18a) 

q^  -  q^  =  -r(x)p(x,t)  .  (18b) 


The  recursions  (17)  -  (18)  constitute  the  fast  Cholesky  recxarsions  [20] , 
and  have  also  been  called  the  downward  continuation  recursions  by  Bube 
and  Burridge  [21] . 

The  initial  data  for  these  recursions  are  the  measured  waves  p(0,t) 
and  q(0,t).  The  algorithm  (17)- (18)  can  be  viewed  as  using  a  layer  stripping 
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principle  to  identify  the  parameters  of  the  scattering  mediiim.  Thus, 
assume  that  the  waves  p(x,t)  and  q(x,t)  at  depth  x  have  been  confuted. 

The  reflectivity  function  r(x)  is  obtained  from  (17)  and  is  used  in  (18) 
to  compute  the  waves  p(x+A,t)  and  q(x+A,t)  at  depth  x+A.  The  effect  of 
the  recursions  (17)  -  (18)  is  therefore  to  identify  and  then  strip  away 
the  layer  [x,x+A). 

The  Schur  Recursions 

The  medixmi  can  also  be  reconstructed  by  using  the  recursions  (3) 
for  the  transformed  waves  p(x;a))  and  q(x,(jo)  with  the  expression 

r(x)  =  2q(x,x) 

=  lim  2jco  exp  jux  q(x,a))  (19) 

where  we  have  assumed  that  the  waves  p  and  q  have  the  form  (16) .  The 
reciirsions  (3)  ,  (19)  constitute  the  frequency  domain  counterpart  of  the 
time-domain  recursions  (17) -(18). 

An  alternate  method  is  to  consider  the  left  reflection  coefficient 

6(x,a3)  =  q(x,w)/p(x,w)  (20) 

which  is  associated  to  the  section  of  the  scattering  medi;mi  extending  over 
[x,°°)  .  The  expression  (20)  assumes  that  the  medium  is  probed  from  the 
left  and  that  no  wave  is  incident  from  the  right.  By  using  the  recursions 
(3)  for  p  and  q,  we  find  that  6(x,a))  satisfies  the  Riccati  equation 

=  2jwR  +  r (x)  (R  -1)  (21) 

and  the  initial  value  theorem  can  be  used  (see  [20])  to  show  that 
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r{x)  =  lim  2jto:6 (x,a)) .  (22) 

When  (22)  is  substituted  in  (21) ,  the  Riccati  equation  (21)  can  be 
propagated  autonomously  and  provides  a  way  of  reconstructing  the 
reflectivity  f\anction  r.  The  initial  condition  for  these  rec-ursions  is 

:&{0,a))  =  6(0))  ,  (23) 

and  is  obtained  from  the  measured  waves  p(0,0))  and  q(0,O))  ,  or  from  direct 
frequency-domain  measvirements  of  the  reflection  coefficient  6(0)). 

It  should  be  emphasized  that  (17)- (18),  (3)  and  (19),  and  (21)- (22) 
are  all  just  different  versions  of  the  same  algorithm  and  are  thus  inter¬ 
changeable.  In  this  paper,  we  will  refer  to  (17)- (18)  as  the  fast  Cholesky 
recursions,  and  to  (3)  and  (19)  or  to  (21)- (22)  as  the  Schur  algorithm, 
but  all  of  these  are  just  different  forms  of  the  Schur  algorithm.  Note 
that  the  Riccati  equation  (21)  for  the  reflection  coefficient  6  is  well- 
known  in  scattering  theory  [22]  and  is  direct  consequence  of  the  rules 
of  composition  of  scattering  layers  [17] .  This  equation  was  in  fact 
used  by  Gjevick  et  al.  [23]  to  develop  an  iterative  method  for  reconstructing 
the  reflectivity  function  r(x).  What  distinguishes  the  Schur  algorithm 
from  these  results  is  the  observation  that  the  relation  (22)  can  be  used 
to  compute  6(x,a))  recursively  for  increasing  values  of  x. 

The  recursions  (21) -(22)  are  the  continuous  version  of  an  algorithm 

obtained  by  Schur  [l]-[2]  for  testing  the  boundedness  of  a  function  R(z) 

which  is  analytic  outside  the  unit  disk.  Given  R(z) ,  Schur  showed  that 

|r(z)I  <  1  outside  the  unit  disk  if  and  only  if  the  reflection  coefficients 

r  obtained  from  the  recursions 
n 
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R  (z)-r 

z(l-r  R  (z)) 
n  n 


Rq(z)  =  R(z) 


(24a) 


r  =  lim  R  (z)  (24b) 

n  ,  n 

zr^oo 

are  such  that  r  <1.  Some  recursions  similar  to  (24)  can  in  fact  be 
n  - 

obtained  by  performing  a  backwards-difference  discretization  of  the 
Riccati  equation  (21).  Similarly  the  fast  Cholesky  recursions  (17)- (18) 
were  used  in  [8]  to  perform  the  causal-anticausal  factorization  of  a 
Toeplitz  operator  and  constitute  the  continuous  counterpart  of  a  discrete 
algorithm  which  was  obtained  in  [5] -[7]  to  construct  the  LDU  factorization 
of  a  Toeplitz  matrix. 

Finally,  it  is  worth  noting  that  the  layer  stripping  principle  which 
was  used  here  to  solve  the  inverse  scattering  problem  for  the  two- 
component  wave  system  (1)  can  also  be  used  for  other  physical  models 
of  scattering  media.  Some  inverse  scattering  techniques  based  on  a  layer 
stripping  principle  were  in  fact  developed  in  [24],  [25]  for  the  telegrapher's 
equation,  and  in  [26]  for  the  Schrodinger  equation  (see  also  [20]).  How¬ 
ever,  instead  of  considering  separately  the  inverse  scattering  problem 
for  a  Schrodinger  equation, we  can  use  the  results  developed  above  for 
two-component  wave  equations. 

The  Schrodinger  Equation 

Consider  the  equation 

y  -Y^4.  =  V(x)y(x,t)  (25) 

XX  tt 

which  is  associated  to  an  elastically  braced  string  [27],  where  y(x,t) 
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denotes  the  displacement  of  the  string  at  point  x  and  time  t,  and  V  (x)  is 
the  elasticity  constant  at  x.  We  assume  that  V(x)  is  localized, 
i.e.  V (x)  =  0  for  x<0  and 

OO 

I  (1+x)  |v(x)  |dx  <  , 

0 

so  that  as  x<0  and  x-><» 


y(x,t)  =  y^(x-t)  i-y^Cx+t) 


(26) 


is  the  superposition  of  rightward  and  leftward  propagating  waves, 
the  Fourier  transform  of  (25)  yields  the  Schrodinger  equation 


y  +  (w  -  V(x))y(x,cj)  =  0  . 

XX 


Taking 


(27) 


Then,  the  inverse  scattering  problem  for  this  equation  is  expressed  in 

terms  of  the  solution  y_  (x,a))  and  $  (x,w)  such  that 

L  K 


and 


-jcox  ,  ,  ncox 

,  e  +  R  (w)e-^ 

Xj 


yj^(x,a))  = 


T  (a})e 
L 


“Dcox 


for  x<0 


as  x-^ 


(28) 


yj^(x,a)) 


!T^(a))e^“^ 


for  x<0 


as 


(29) 


which  define  the  scattering  matrix  S(ui)  associated  to  V  (x)  .  These  solutions 
correspond  to  the  case  when  the  string  is  probed  from  the  left  or  from  the 
right  by  an  impulsive  wave.  The  problem  is  to  determine  V(x)  given  the 


-15- 


reflection  coefficient  function  or 

Several  solutions  of  this  prcblem  based  on  the  Gelfand-Levitan 
procedure  [18]- [19],  [28]  or  on  trace  formulas  [29]  have  been  proposed. 

We  will  now  show  that  the  Schur  algorithm  provides  a  solution  which  is 
computationally  quicker.  To  see  this,  note  that  by  taking  the  derivative 
of  the  two-conponent  system  (3)  with  respect  to  x,  we  obtain  the  matrix 
Schrodiiiger  equation 


liS  i- 

"2 
r  -  r 

X 

\ 

p(x,0)) 

2 

-r  r 

q  (x  ,(j0) 

\  ^ 

X 

/ 

where  denotes  the  2x2  identity  matrix.  By  making  the  change  of  variable 


(x,0)) 

=  p(x,0))  +  q(x,w) 

Ola) 

y2  (x/W) 

=  p(x,0))  -  q(x,a)) 

(31b) 

this  equation  can  be  decoupled  into  two  scalar  Schrodinger  equation 


y^^^  +  (x) )y^ (x,w)  =  0 

^2xx  ^  (x)  )y2  (x,a))  =  0 

where 

V^(x)  =  r^(x)  -  r^(x) 

V  (x)  =  r^(x)  +  r  (x) 


(32a) 

(32b) 


(33a) 

(33b) 


In  addition,  we  observe  from  (31)  and  from  the  definition  of  the  scattering 
matrix  3(0))  of  the  two-component  system  (3)  that  the  scattering  matrix 
associated  to  V^(x)  is  identical  to  that  of  (3),  and  that  the  scattering 
matrix  82(0))  associated  to  given  by 
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(W)  = 


f  (0))  -  (OJ) 
Xi  K 

-6^(0))  T^(U)) 


ZS(W)E  , 


(34) 


i.e.  it  is  obtained  by  changing  the  sign  of  the  reflection  coefficients 
and  of  (3). 

Consequently,  given  a  potential  V (x) ,  we  can  always  view  its  left 


reflection  coefficient  R^.  (w)  as  arising  from  a  two-couponent  system  such 

Ij 

as  (3) .  Then,  given  (w)  or  the  impulse  response  R  (t) ,  we  can  use 

the  Schur  or  fast  Cholesky  recursions  to  reconstruct  the  reflectivity  function 

r(x),  which  in  turn  can  be  used  to  recover  V(x)  from  the  relation  (33a). 

The  relation  (33a)  is  known  in  soliton  theory  as  the  Miura  transformation 
[10] ,  [27]  ,  and  it  maps  solutions  of  the  modified  Korteweg-de  Vries 
equation  into  solutions  of  the  Korteweg-de  Vries  equation. 

If  the  potential  V(*)  extends  only  over  a  finite  interval  [0,L], 
the  interval  [0,L]  may  be  divided  into  N  sxjbintervals  of  length  A  =  L/N, 
and  the  Schur  and  fast  Cholesky  recursions  may  be  discretized  accordingly. 

It  is  shown  in  [20]  ,  [21]  that  the  resulting  procedure  requires  only 

2  3 

0  (N  )  operations  to  recover  r(*)  and  V(-),  instead  of  0 (N  )  if  we  dis¬ 
cretize  the  Gelf and- Levitan  equation  and  solve  the  resulting  system 
of  linear  equations.  The  Schur  and  fast  Cholesky  algorithms  are  therefore 
quite  efficient.  Note  however  that  if  we  exploit  the  structtire  of  the 
Gelfand-Levitan  equation  and  use  the  Levinson  recursions  to  solve  this 
equation  [20] ,  we  obtain  a  reconstruction  procedure  which  is  as  efficient 


as  the  Sch'ur  algorithm. 
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lll*  The  Lossless  'Nonrtuniform  Transmission  Line 


In  this  section  we  study  the  inverse  problem  for  the  lossless  non- 
uniform  transmission  line,  and  show  that  its  solution  is  given  by  the  Schur 
or  fast  Cholesky  algorithms  (see  [30]  for  an  earlier  solution  of  this  problem) . 
In  the  process,  we  give  a  scattering  interpretation  of  transmission  line 
phenomena  such  as  waves,  reflections,  impedances. 

Consider  an  infinitesimal  section  of  length  A  of  a  lossless  non-uniform 
transmission  line.  Such  a  section  is  illustrated  in  Fig.  3.  Note  that  L(x) 
and  C(x)  represent  inductance  and  capacitance  per  unit  length,  i.e.  they  are 
distributed  quantities.  Writing  equations  for  Fig.  3,  we  have 
v(,x,t)  =  Li^A  +  v(x  +  A,  t) 


i(x,t)  =  Cv^A  +  i(x  +  A,  t) 


(35) 


Dividing  by  A  and  letting  A  0,  we  obtain  the  telegrapher's  equations 
V  +  L(x)i^  =  0 

X  t 

i  +  C(x)v^  =  0 
X  t 


(36) 


which  also  arise  in  acoustics  [24]  and  in  studies  of  the  human  vocal  tract 

[25],  [31]  under  the  assumption  of  losslessness. 

For  a  uniform  line,  it  is  well  known  (see  [32])  that  (36)  admits  wave 

solutions,  and  that  for  such  waves  the  ratio  of  the  amplitudes  of  the  voltage 

1/2 


and  current  is  the  characteristic  impedance  =  (L/C) 


Since  the  quantities 
p  and  q  appearing  in  the  two- component  wave  equations  must  be  dimensionally 
equivalent,  this  suggests  defining  for  the  non-uniform  line  the  dimensionally 
equivalent  variables 


-1/2 

V(x,t)  =  Zq  v(x,t) 
I(x,t)  =  Zq'^^  i(x,t) 


(37) 
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1/2 

where  Z^Cx)  =  (Ii(x)/C(x))  .  Substituting  (37)  in  (36)  yields 


(38) 


In  order  to  make  the  dependent  variables  x  and  t  dimensionally  equivalent, 
we  replace  x  with  the  travel  time  z  defined  by 


z  (x) 


(L(u)  C(u))^^^  du 


(39) 


0 

1/2  , 


Since  (L(x)  C (x) )  “  is  the  local  wave  speed  at  x,  z (x)  is  the  time  required 

for  a  wave,  starting  at  x  =  0,  to  reach  position  x.  Making  the  additional 
change  of  variables 


p(z,t)  =  —  (V(z,t)  +  I(z,t)) 


q(z,t)  =  —  (V(z,t)  -  I(z,t))  , 


(40a) 

(40b) 


and  defining  the  reflectivity  function 


(41) 


we  obtain  the  two-component  wave  equations  (1) .  The  relations  (40)  provide 
an  interpretation  of  the  right  and  left  propagating  waves  in  terms  of  the 
normalized  voltage  and  current. 

Interpretation  of  the  Reflection  Coefficient 


Suppose  a  xmiform  transmission  line  is  terminated  with  a  load  Z  .  Then  a 

Xj 

/\ 

wave  travelling  down  the  line  will  be  reflected  back  by  the  load.  Define  R(U))  , 
the  reflection  coefficient  for  the  load,  to  be  the  ratio  of  the  Fourier  transforms 
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of  the  primary  and  reflected  voltage  waves ,  at  the  frequency  w.  It  is  easy 
to  show  (see  132] )  that 


R(w) 


V  (w) 
KEFL 

V  (W) 
PRIM 


(42) 


For  the  non-uniform  transmission  line  considered  here,  since  there  is  a 
one-to-one  correspondence  between  position  x  and  travel  time  z,  we  will  use  x 
instead  of  z  in  the  qualitative  analysis  to  follow.  Then,  at  point  x  on  the 
line,  the  load  perceived  due  to  all  of  the  line  to  the  right  of  x  is  (see  Fig.  4) 


z  (,x,a))  =  v(x,aj) /i (x,a))  .  (43) 

Xj 

By  substituting  this  expression  in  (42) ,  we  find  that  for  the  non-uniform 
transmission  line,  the  reflection  coefficient  at  point  x  is 


Vi  -  Z  (x) 

R(x,lj^l  =  - 

v/i  +  Zq(x) 


/S 

V/I  -  1 
V/I  +  1 


=  q(x,co)/p  (x,a)) 


(44) 


This  is  precisely  the  expression  (20)  for  the  left  reflection  coefficient  of  the 
section  of  the  two-component  system  (1)  extending  over  [x,'»)  . 

We  see  therefore  the  meaning  of  R(z,co)  .  For  a  given  point  x  on  the  line, 
and  any  given  frequency  Ci),  it  is  the  ratio  of  the  reflected  and  primary  voltage 
waves,  with  the  reflection  due  to  the  inhomogeneity  of  the  line  at  x.  Prom 
Section  II,  we  know  that  R(z/a))  satisfies  the  Riccati  equation  (21)  ,  and  that 

/N 

r(x)  may  be  forad  from  R(x,a))  by  using  (22)  .  Also  note  that  if  the  line  is 

dZo 

locally  uniform  at  point  x.,  q —  (x„)  =  0  hence  r(x,,)  =  0  and  no  reflection 

0  dx  0  0 
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occurs .  Reflections  occur  only  where  the  line  is  inhomogeneous , 


Inverse  Problem 


Suppose  now  that  the  line  characteristics  L(x)  and  C(x)  are  unknown  and 
that  we  want  to  determine  them  from  the  measured  impedance  Z  (u)  =  Z  (0,0)).  This 

Ij 

problem  arises  not  only  when  we  want  to  find  the  characteristics  of  an  existing 
transmission  line,  but  also  if  we  want  to  synthesize  a  transmission  line  with 
prescribed  impedance  Z  (o))  ,  It  is  assumed  here  that  we  have  access  to  only  one 
end  of  the  line.  The  line  characteristics  can  be  partially  reconstructed  as 
follows.  First,  set  2^(0)  =  1  and  consider  the  reflection  coefficient 


r'(w) 


Z(C0)  -  1 
Z{W)  +  1 


(45) 


Then,  run  the  Schur  recursions  (21)  -  (22)  ,  using  R(W)  as  initial  condition, 
to  obtain  r(z).  Alternately,  we  may  compute  the  inverse  Fourier  transform  R(t) 
of  R(w),  and  use  the  fast  Cholesky  recursions  (17)  -  (18)  to  obtain  r(z) .  Given 
r(z),  the  expression 


Z  (z) 
0 


exp  2 


•  z 

r(u)  du 

■'  0 


(46) 


1/2 

enables  us  to  recover  the  characteristic  impedance  Zq(z)  =  (L(z)/C(z))  as  a 
function  of  the  travel  time  z.  However,  we  cannot  reconstruct  L(x)  and  C(x) 
separately  as  functions  of  the  position  x. 

The  same  difficulty  will  appear  in  Section  IV  for  the  inverse  seismic  problem 
of  geophysics,  except  that  in  this  case  we  will  be  able  to  use  an  additional  degree 
of  freedom,  the  angle  of  incidence  of  the  probing  waves,  in  order  to  reconstruct 


the  medium  completely. 
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IV.  The  Inverse  Seismic  Problem 

In  this  section  we  examine  and  solve,  using  the  Schur  algorithm,  the 
inverse  problem  for  a  one-dimensional  acoustic  medium  probed  with  plane  waves. 

We  will  consider  first  the  case  of  plane  waves  at  normal  incidence,  and  then  the 
offset  problem  in  which  the  probing  waves  come  in  at  an  angle ,  as  shown  in  Pig .  5 . 

For  the  normal  incidence  case,  the  equations  we  will  obtain  are  almost 
identical  to  those  of  last  section,  and  even  though  the  physical  situations  are 
quite  different,  the  inverse  problem  is  the  same  as  for  the  lossless  non-uniform 
transmission  line. 

A  simple  transformation  will  allow  the  offset  problem  to  be  solved  by  using 
the  same  procedure  as  for  the  normal  incidence  problem.  By  probing  the  medium 
at  two  different  angles,  it  will  be  shown  that  the  medium  density  and  velocity 
profiles  can  be  reconstructed  separately  as  functions  of  depth.  The  use  of  the 
Schur  algorithm  to  solve  the  offset  problem  has  not  to  our  knowledge  appeared  in 
the  literature. 

The  Normal  Incidence  Problem 

The  problem  to  be  considered  in  this  section  corresponds  to  the  case  when 
the  angle  of  incidence  0  =  0  in  Fig.  5,  The  acoustic  medium  that  we  want  to  re¬ 
construct  is  constituted  of  a  homogeneous  half-space  with  known  density  and 
sound  speed  eg  extending  over  x  <  0,  and  of  an  inhomogeneous  half-space  with 
unknown  density  p(x)  and  unknown  local  sound  speed  c(x)  extending  over  x  >  0. 

For  convenience,  we  assume  that  inhomogeneities  are  localized,  so  that  for  x  >  L 
the  density  and  sound  speed  c^  are  constant.  Physically,  the  region  x  <  0 
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corresponds  to  the  air  or  ocean  located  above  the  mediiim  to  be  probed,  and 
the  region  x  >  L  corresponds  to  the  substrate  or  bedrock  located  below  it.  An 
impulsive  plane  pressure  wave,  propagating  downwards,  is  incident  upon  the  in¬ 
homogeneous  region  at  t  =  0,  and  the  reflections  or  reverberations  making  their 
way  back  to  the  surface  of  the  inhomogeneous  medium  are  measured  either  at  the 
surface  of  this  medium  (land  case)  or  above  the  surface  (marine  case) .  Our  goal 
is  to  botain  profiles  of  p (x)  and  c (x)  as  functions  of  depth.  The  presentation 
will  follow  that  of  [33]  -  [36] . 

The  two  basic  equations  we  start  with  are  the  acoustic  equation  [33] ,  [35] 


(47) 


(48) 


where  w(x,t)  and  P(x,t)  denote  respectively  the  displacement  and  pressure  (negative 
stress)  at  depth  x. 

The  first  step  is  to  change  variables  from  depth  x  to  travel-time  z (x) ,  which 
is  the  time  it  takes  for  a  wave  starting  at  the  surface  of  the  inhomogeneous  medi-um 
to  reach  depth  x.  (Recall  a  similar  definition  in  the  last  section) .  Thus  we 
have 


z  (x)  = 


' 

du/c (u) 

0 


(49) 


By  substituting  (49)  inside  (47)  -  (48) ,  and  defining 


Z(z)  =  p(z)c(z)  =  characteristic  impedance  (50a) 

v(z,t)  =  w^(z,t)  =  particle  velocity 


(50b) 
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we  obtain  the  system 

P  =  -Z(z)v, 
z  t 

V  =  -  Z~^(z)P^ 

z  t 

which  in  the  Fourier  domain  takes  the  form 

/V  /s 

p  =  -  j  coz  ( z  )  V 
z 

V  =  -  j  Z  (z)P 
z 

Then,  if 

-1/2 

'i'(z,t)  =  Z  (z)  P(z,t)  =  normalized  pressure 
1/2 

(j)(z,t)  =  Z  (z)  v(z,t)  =  normalized  velocity 


(51) 

(52) 

(53a) 

(53b) 


and  if  we  make  the  change  of  variables 

p(z,t)  =  I  {'?  +  (f>) 

q(z,t)  =  -I  ('i'  -  4)) 


(54) 


the  system  (51)  can  be  transformed  into  the  two- component  wave  system 


p  +  p  =-r{z)q(z,t) 
z  t 

q^  -  q^  =  -  r(z)p(z,t) 

where  the  reflectivity  function  r(z)  is  given  by 
r(z)  =  —  ^  Z(z) 


(55) 


(56) 


The  definition  of  the  normalized  variables  ¥(z,t)  and  4)(z,t)  and  of  the 
waves  p(z,t)  and  q(z,t)  is  identical  to  the  one  we  used  in  last  section  for  the 
normalized  voltage  and  current  V(z,t)  and  I(z,t)  and  the  associated  right  and 
left  going  waves.  The  inversion  problem  for  the  1-D  acoustic  medium  probed  by 
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plane  waves  at  normal  incidence  is  therefore  the  same  as  that  for  the  nonuniform 
lossless  transmission  line. 

The  data  which  is  used  in  the  inversion  is  obtained  by  sending  a  downward 
impulsive  plane  pressure  wave  which  is  normally  incident  upon  the  inhomogeneous 
half-space  at  t  =  0.  Then,  for  the  land  case,  we  measure  the  particle  velocity 
v{0,t)  at  the  earth's  surface  as  a  function  of  time.  Since  the  difference  in 
density  between  the  air  and  the  earth  is  very  large,  the  earth's  surface  acts 
like  a  free  surface,  i.e.  the  pressure  on  it  is  zero  for  positive  times.  We 
can  therefore  express  the  pressure  and  velocity  on  the  surface  as 


P(0,t)  =  PQ'5(t)  (57a) 

v(0,t)  =  VQ(6(t)  +  2fi(t)u(t))  (57b) 

where  Pq/'Vq  =  Z(0).  After  normalization,  the  downgoing  and  upgoing  waves  p(0,t) 
and  q(0,t)  take  the  form 


p(0,t)  =  6(t)  +  h(t)u(t) 
q(0,t)  =  -  h(t)u(t)  , 


(58) 


and  the  fast  Cholesky  recursions  (17)  -  (18)  or  the  Schur  recursions  (3)  and  (22) 
can  be  applied  to  these  waves  to  reconstruct  the  reflectivity  function  r(z). 

For  the  case  of  a  marine  seismogram,  the  reflected  pressure  wave  R(*)  is 
measured  at  some  point  inside  the  homogeneous  half-space  x  <  0.  The  pressure  and 
velocity  in  this  half-space  are 


P(x,t)  =  P^(5(t-x/cQ)  +  R(t  +  x/Cq)  u(t  +  x/Cq)) 
v(x,t)  =  VQ(6(t-x/cQ)  -  R(t  +  x/Cq)  u(t  +  x/Cq) ) 


(59) 
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so  that  at  the  surface  of  the  inhomogeneous  medium  {the  ocean  floor) ,  the 
downgoing  and  upgoing  waves  are  given  by 

p(0,t)  =  6 (t) ,  q(0,t)  =  R(t)u(t)  .  (60) 

These  can  then  be  used  to  reconstruct  the  reflectivity  function  r(z). 

Given  r(z),  the  impedance  Z(z)  can  be  obtained  by  using  (46).  Thus, 
as  in  the  case  of  the  nonuniform  lossless  transmission  line,  the  best  we  can 
do  is  to  reconstruct  Z(z)  =  pc  as  a  function  of  the  travel  time  z.  We  cannot 
recover  p (x)  and  c(x)  separately  as  functions  of  depth.  The  procedure  consisting 
of  using  the  Schur  recursions  to  reconstruct  the  impedance  Z(z)  is  commonly 
called  dynamic  deconvolution.  It  was  developed  first,  using  the  discrete 
Schur  algorithm,  for  the  case  of  a  layered  medium  divided  up  into  homogeneous 
layers  of  equal  travel-time  (the  so-called  Goupillaud  model)  and  it  is  described 
in  [9]. 

Another  feature  of  this  reconstruction  method  is  that  the  quantities  p(x,t) 
and  q(x,t)  represent  respectively  downgoing  and  upgoing  waves.  So  when  we  run 
the  fast  Cholesky  or  Schur  recursions  on  the  experimental  data,  we  are  decomposing 
the  pressure  and  particle  velocity  at  each  depth  into  a  superposition  of  upgoing 
and  downgoing  waves.  Thus,  we  gain  not  only  information  about  the  medium 
parameters,  but  also  information  about  what  is  happening  to  the  medium.  This 
could  prove  useful  in  evaluating  how  realistic  the  model  is. 

The  Offset  Problem 

We  now  consider  the  problem  in  which  an  impulsive  plane  pressure  wave  is 
obliquely  incident  on  the  medium  at  an  angle  0  to  the  vertical,  as  shown  in 
Fig.  5.  In  this  case,  the  impulse  response  R(t,  y;  0)  is  a  function  of  the 
horizontal  coordinate  y  (in  the  normal  incidence  case  there  is  of  course  no 
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horizontal  variation) .  We  will  once  again  obtain  a  dynamic  deconvolution 
procedure  that  uses  the  Schur  algorithm  to  recover  an  impedance  as  a  function 
of  travel  time,  although  this  impedance  differs  slightly  from  (50a) .  However, 
running  the  offset  experiment  twice  for  two  different  angles  of  incidence 
0  =  6^,  ©2  will  allow  us  to  recover  p (x)  and  c{x)  separately,  as  functions 
of  depth,  a  significant  improvement  over  the  normal  incidence  experiment. 

The  problem  is  set-up  as  in  [37]  -  [39] . 

An  impulsive  plane  pressure  wave  6{t-(x  cos  0  +  y  sin  Q)/c^)  is 
incident  at  an  angle  0:  from  the  vertical,  where  x  and  y  are  as  in  Fig.  5. 

The  Fourier  transform'  of  this  wave  is  P^e  x  y^  ,  where  ^  ^0 

and  k  =0)  sin  6/c„  are  the  vertical  and  lateral  wavenumbers  in  the  upper 

y  0 

homogeneous  half-space.  For  the  case  of  a  marine  seismogram,  the  pressure 
field  for  x  <  0  is  therefore 

^  -jk  y  -jk  X  ^  jk  X 

P(x,y,a);0)  =  P^e  ^  (e  ^  +  R((jJ;0)e  ^  (61) 

(compare  this  to  the  Schrodinger  equation  boundary  condition  (28)).  This 
shows  that  in  the  time  domain  the  impulse  response  R(t,y;0)  has  the  form 

R(t,y;0)  =  R(t-y  sin  ©/c^;©)  (62) 

where  R(t;9)  is  the  inverse  Fourier  transform  of  the  reflection  coefficient 
R(a3;6).  This  form  is  also  valid  for  the  case  of  a  land  seismogram.  Thus, 
in  theory  it  should  only  be  necessary  to  measure  R(t,y;0)  at  a  single 
surface  point  (e.g.  y  =  0) .  However,  in  practice  we  need  to  take  data  for 
a  range  of  y  and  filter  or  stack  it  to  the  form  (62) .  This  is  because  any 
real-world  impulsive  wave  can  only  be  locally  planar,  while  the  form  (61)  of 
the  pressure  field  assumes  an  incident  plane  wave  of  infinite  extent. 
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With  y  dependence  added,  the  acoustic  and  stress-strain  equations  (47) , 
(48)  become 


=  -  Vp(x,y,t)  (63a) 

2 

P(x,y,t)  =  -  pc  V  •  w(x,y,t)  (63b) 


where  the  displacement  w(x,y,t)  is  now  a  vector.  Then,  if  v  =  w^  is  the 

X  V 

particle  velocity,  and  if  v  and  are  its  vertical  and  lateral  components, 
by  taking  Fourier  transforms  we  obtain 


P^  =  -  jcop(x)v 

^  /sV 

P^  =  -  3a}p(x)v 

,  ,  2 ,  .  ,^x  ,  ^y,  .  r; 

p(x)c  (x)  (v  +  V  )  =  -  3CjOP 
X  y 


(64a) 

(64b) 

(64c) 


Since  the  medi\im  properties  vary  only  with  depth  x,  the  horizontal 


wavenumber  is  preserved,  and  we  may  write,  as  in  [37] 


“jk.y 


P(x,y,a))  =  'rT(x,co)e 


(65) 


S\abstituting  this  expression  in  (64b)  yields 

v^  =  (sin  6/p(x)Cq)P  ,  (66) 

and  by  inserting  (66)  in  (64c)  and  using  Snell’s  law 

sin  0(x)/c(x)  =  sin  6/Cq  =  ray  parameter  (constant)  (67) 

where  6 (x)  is  the  local  angle  that  a  ray  path  makes  with  the  vertical, 
we  get 

p(x)c^(x)v^  =  -  jto  cos^  0(x)P  . 


(68) 
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Now,  define 

c' (x)  =  c (x) /cos  0 (x)  =  local  vertical  wave  speed  (69) 

z(x)  =  du/c' (u)  =  vertical  travel  time  to  depth  x  (70) 

■*0 

z (z)  =  p(z)c'(z)  =  effective  impedance  .  (71) 


Using  (69)  =  (71)  in  (64a)  and  (68)  gives 


P  =  “  jco  Z  (z)  V 
z 


/sv  —  1 

V  =  -  jW  Z  (z)P  , 


and  once  again  defining  the  downgoing  and  upgoing  waves  as 


p(z,y,a))  =  ^  (2  ^'^^P"(z,y,(i))  +  Z^'^^V^(z,y,(iO)) 
q(z,y,a))  =  y(2' (z,y,w)  -  z^^^’(^(z,y,0J) )  , 


(72a) 

(72b) 


(73a) 

(73b) 


we  obtain  the  two-component  wave  system 

=  -  jujp  -  r(z)q  (74a) 

q^  =  -  r(z)p  +  jwq  (74b) 

where  the  reflectivity  function  r(z)  is  given  by  (56).  Here  y  is  fixed  at 

whatever  value  of  y  we  measure  R(t,y;0)  (e.g.  y  =  0) ,  and  6  is  a  parameter  on 

which  all  quantities  depend. 

Note  that  once  again  the  quantities  in  the  wave  system  (74)  are  the 
Fourier  transform  of  the  downgoing  and  upgoing  waves,  so  that  the  vertical 
motion  of  the  medium  is  again  decomposed  into  upgoing  and  downgoing  waves. 
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Furthermore,  for  both  the  case  of  a  land  and  marine  seismogram,  these 
waves  have  the  form  (16) ,  so  that  we  can  use  the  fast  Cholesky  or  Schur 
recursions  to  reconstruct  the  effective  impedance  Z(z)  given  by  (71). 

Now,  suppose  that  the  offset  experiment  is  run  twice,  for  two  angles 
of  incidence  6=0,  0„.  Two  different  imepdances  Z  (z  )  and  Z  (z  )  are 

•i-  ^  X  X  ^ 

obtained,  which  are  functions  of  two  different  travel  times  z^  and  z^-  From 
(69)  -  (71) ,  we  get 


cos  9  (x)/cos  0  (x)  =  Z  (z  )/Z  (z  ) 

X  2.  2.  2  L 


(75) 


which  can  be  solved  or  integrated  numerically  to  obtain  the  monotone  increasing 
function  z^  =  ^^(z^).  This  unables  us  to  express  Z^  and  Z^  as  functions  of 
the  same  travel  time  z^.  Then,  using  (67)  and  (69)  -  (71) ,  we  can  reconstruct 
separately  as  functions  of  the  known  impedances  Z^  and  Z^. 
Finally,  inverting  (70)  and  using  (69)  gives  the  effective  travel  time  z^ (x) , 
yielding  p(x)  and  c(,x)  separately,  as  functions  of  depth. 

This  reconstruction  procedure  has  only  been  sketched,  since  it  has  nothing 
to  do  with  the  Schur  algorithm;  it  is  presented  in  more  detail  in  [38].  Note 
however  that  the  dynamic  deconvolution  method  described  above  to  compute 
the  impedances  2^(z^)  i  =  1,  2  is  new.  An  alternate  reconstruction  procedure 
was  described  in  [39]  which  recovers  p(x)  and  c(x)  recursively  by  operating 
directly  on  the  need  to  compute  the  impedances  Z^(z^)  i  =  1,  2  as  a 
preliminary  step. 

The  reason  that  the  profiles  p Cx)  and  c(x)  recursively  by  operating  directly 
on  the  waves  associated  to  0^  and  6^ ,  thereby  obviating  the  need  to  compute 
the  impedances  Z^(z^)  i  =  1,  2  as  a  preliminary  step. 
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The  reason  that  the  profiles  p(x)  and  c(x)  can  be  recovered  separately 
for  the  oblique  incidence  problem,  but  not  for  the  normal  incidence  problem, 
is  that  by  running  the  oblique  experiment  twice  information  has  been  gained 
along  two  different  ray  paths.  This  option  is  not  available  for  the  normal 
incidence  problem. 

Note  that  along  any  given  ray  path  the  ray  parameter  (67)  is  constant, 
so  that  unless  the  angle  of  incidence  0  is  less  than  the  Critical  angle 
sin  ^  (c^/max  c(x)),  the  angle  0(x)  will  become  imaginary  at  some  depth. 
Physically,  this  situation  results  in  evanescent  waves,  in  which  the  pressure 
field  decays  exponentially  with  depth.  This  causes  no  difficulty  in  the 
Schur  algorithm  until  the  ray  path  becomes  horizontal,  prior  to  turning 
back  up.  When  this  turning  point  is  reached,  |r(x)  |  «>.  However,  since 

no  new  information  can  be  gained  beyond  the  turning  point,  the  exceeds  a 
pre-set  value.  For  the  reconstruction  procedure  described  above,  this  means 
that  we  can  recover  p(x)  and  c (x)  only  until  a  turning  point  occurs  in  either 
of  the  two  oblique  incidence  experiments. 


V.  Linear  Estimation  of  a  Stationary  Stochastic  Process 

In  this  section  the  problem  of  finding  the  linear  least-squares  estimate 
of  a  stationary  stochastic  process  given  some  observations  of  this  process  over 
a  finite  interval  is  posed  as  an  inverse  scattering  problem,  and  solved  using 
the  Schur  algorithm.  This  formulation  of  the  estimation  problem  for  a  stationary 
stochastic  process  is  due  to  Dewilde  and  his  coworkers  [4] ,  [15] ,  [16] ,  [40]  . 


The  basic  problem  to  be  considered  is  as  follows.  Let 
y(.t)  =  z{t).  +  v(t) 


(76) 
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be  some  observations  of  a  zero-mean  stationary  stochastic  process  z(-) 
with  covariance 

E[z(t)z(.s)]  =k(|t-sl),  {77a) 

where  v(-)  is  a  white  noise  process  with  unit  intensity,  i.e. 

E[v(t)v(s)]  =  6(t-s)  .  (77b) 


We  assume  that  z(-)  and  v(-)  are  uncorrelated  and  that  k(*)  e  L^[0,°°),  so 
that  its  Fourier  transform 


k(0J)  = 


k(t)  exp- jut  dt 


(78) 


exists.  In  this  case,  the  spectral  density  of  y(*)  is  w(a))  =  1  +  k(a)) 
+  k(-a))  . 

Given  the  Hilbert  space 


Y (t;  x)  =  H (y (t+s) ,  -x  <  s  <  x) 


(79) 


spanned  by  the  observations  over  the  interval  [t-x,  t+x] ,  our  objective  is 
to  compute  the  forwards  and  backwards  linear  least-squares  estimates  of  z 
at  the  endpoints  of  this  interval.  These  estimates  can  be  denoted  as 


z (t+x)  Y(t;  x) ) 

z(t-xlY{t;  x) ) 


■X 

A(x;  u)  y(t+u)du 

-X 

B(x;  u)  y(t+u)du  , 

'  -X 


(80a) 


(80b) 


where  A(x;  •)  and  B(x;  •)  are  the  optimal  forwards  and  backwards  prediction 
filters,  respectively.  Note  that  since  the  process  z(-)  is  stationary  the 
filters  A(x;  •)  and  B{x;  •)  do  not  depend  on  t,  the  center  of  the  interval 
[t-x,  t+x] .  Then,  if  the  forwards  and  backwards  residuals  are  defined  as 
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e(t,x)  =  y(t+x)  -  z(t+x|Y(t;  x) )  (81a) 

b(t,x)  =  y(t-x)  -  z(t-x|Y{t;  x) )  ,  (81b) 

by  using  the  orthogonality  property 

e(t,x),  b(t,x)  J_Y(t;  x)  (82) 

of  linear  least-squares  estimates,  we  find  that  the  filters  A(x;  •)  and 
B(x;  •)  satisfy  the  Wiener-Hopf  equations 


A(x; 

s) 

+ 

• 

•X 

A(x; 

-x 

u) 

k ( 1 u-s 1 ) 

du 

=  k(x-s) 

(83a) 

B(x; 

s) 

+ 

B  (x; 

-X 

u) 

k ( 1 u-s 1 ) 

du 

=  k  (x+s) 

(83b) 

with  -X  <  s  <  X. 

3  3  3  3 

Applying  the  operators  ^  and  to  (83a)  and  (83b)  respectively, 

and  using  the  linearity  of  the  resulting  equations  yields  the  Krein-Levinson 
recursions  [15] ,  [41] 

8  9 

^■9x  ^  (84a) 

(v - S-)  B(x;  s)  =  -  r(s)  A(x;  s)  (84b) 

ox  os 

with  -  X  <  s  <  X,  and  where 

r (x)  =  2A(x;  -x)  =  2B(x;  x)  (85) 

is  the  reflectivity  function.  The  last  identity  in  (85)  is  obtained  by  noting 
from  a  time-reversal  arg'ument  that  B(x;  s)  B  A(x;  -s)  .  The  Krein-Levinson 
recursions  (84)  have  the  same  form  as  the  fast  Cholesky  recursions.  However, 
as  noted  in  [20] ,  these  two  sets  of  recursions  differ  by  the  fact  that  the 


-33- 


Krein-Levinson  recursions  correspond  to  a  boundary  value  problem  where 
r (x)  is  computed  at  every  step  from 


r(x)  =  2(k(2x) 


A (x;  u)  k (x+u) du) , 


(86) 


whereas  the  fast  Cholesky  recursions  give  rise  to  an  initial  value  problem. 

The  recursions  (84) ,  (86)  can  be  used  to  compute  efficiently  the  filters 

3  “  9 

A(x;  •)-  and  B(x;  ')  •  Furthermore,  if  we  apply  the  operators  7—  +  to  the 

ox  dt 

definition  (81)  of  the  forwards  and  backwards  residuals  e(t,x)  and  b(t,x) 
and  use  the  Krein- Levinson  recursions  (84) ,  we  obtain 


3  3 

(■^  -  •^)  e(t,x)  =  -  r(x)  b(t,x) 

3  3 

(.TT-  +  1^)  b(t,x)  =  -  r(x)  e(t,x) 


(87a) 

(87b) 


This  shows  that  the  residuals  satisfy  a  two-component  wave  system,  where 
e(t,x)  and  b(t,x)  propagate  respectively  leftward  and  rightward,  and  where 
the  waves  at  x  =  0  are  given  by 

e(t,  0)  =  b(t,  0)  =  y(t)  .  (88) 

As  a  consequence  of  this  observation,  the  process  y(t)  can  be  viewed  as  the 
output  of  a  modeling  filter  driven  by  e(t,x)  as  shown  in  Fig.  6a.  This  modeling 
filter  is  obtained  by  aggregating  infinitesimal  ladder  sections  of  the  type 
described  in  Fig.  6.b- 

The  scattering  matrix  associated  to  the  two-component  wave  system  (87) 

can  be  identified  by  noting  that  as  x 

e(t,x)  =  V  (t+x) ,  b{t,x)  =  V  (t-x)  (89) 

F  B 

where  V  ( - )  and  V  ( • )  denote  respectively  the  forwards  and  backwards 
I"  B 

innovations  processes  associated  to  y(*)  [41],  The  processes  V  (•)  and  V  (•) 

T  B 

are  white  noise  processes  and  are  related  to  the  observations  y(’)  through 
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the  identities 


£ 

II 

£' 

<;> 

(90a) 

y(03)  =  F(-03)  V  (03) 

B 

(90b) 

/N  /N  /S 

where  y(w)/  V  t(jO)  and  V  (oj)  denote  formally  the  Fourier  transforms  of 
F  B 

y  { • )  /  V  ( • )  and  V  ( • )  ,  and  where  the  shaping  filter  F  (co)  is  the  outer 
r  B 

spectral  factor  of  w(a3)  ,  i.e. 

w(ta)  =  |f(w)|^  (91) 

on  the  real  axis,  and  F(a))  and  F  ^(co)  are  analytic  in  the  lower  half-plane. 

The  relations  (88)  and  (89)  imply  that  the  scattering  matrix  S  (co) 
satisfies 


/N 

/\ 

Vsto) 

=  S(03) 

y  (03) 

y  (03) 

(  03) 

F 

(92) 


and  by  sxibstituting  (90)  inside  this  relation,  we  obtain  the  identity 


—  • 

“ 

- 

-1 

/N 

F  (-03) 

T 

R„ 

1 

L 

R 

— 

/\ 

“1  ,  V 

1 

R 

T 

F  (03) 

L 

R  _ 

(93) 


for  the  entries  of  S(W).  By  using  the  properties  (11), 
matrix,  this  gives  after  some  algebra 


it  (03) 
Xj 


k'(03) 

1  +  k  (03) 


(12)  of  the  scattering 


(94a) 


(03)  =  T^{03) 


F  (03) 

1  +  k  (03) 


(94b) 


(94c) 


1  +  k  (co) 


F(a3) 

F(-w) 


A,  , 

where  the  left  reflection  coefficient  R  (to)  depends  only  on  the  covariance 

Xj 

A 

data  given  by  k  (to)  . 

Then,  we  observe  that  the  Krein-Levinson  recursions  (84)  and  the  system 
(87)  for  the  residuals  e(t,x)  and  b(t,x)  are  parametrized  entirely  by  the 
reflectivity  function  r(x) .  Consequently,  the  linear  least-squares  estimation 
problem  over  an  arbitrary  finite  interval  will  be  solved  completely  once  we 
reconstruct  rC-).  This  problem  can  be  formulated  as  an  inverse  scattering 
problem  where  R  (co)  is  given  in  the  form  (94a)  ,  and  where  we  want  to  recover 

li 

r  (x)  . 


To  do  so,  one  method  is  to  apply  the  Schur  or  fast  Cholesky  recursions 

directly  to  R  (W)  or  its  inverse  Fourier  transform  R  (t) .  However,  the 
L  Xj 

special  form  of  (94a)  can  be  exploited  by  selecting 


p(0,t)  =  6  (t)  +  k(t)u(t) 
q(0,  t)  =  k  (t)u(t) 


(95) 


as  probing  waves  (see  [8],  [40]),  to  which  we  can  then  apply  the  fast  Cholesky 
recursions.  In  this  case,  as  a  byproduct  of  the  fast  Cholesky  algorithm,  we 
obtain  a  factorization  of  the  covariance  operator  w(t-s)  =  6(t-s)  +  k(|t-s|)  in 
terms  of  causal  times  anticausal  Volterra  operators  [8] .  Furthermore  by  noting 


that 


F(W) 

1  +  k(a)) 

=  S((jO) 

k(a)) 

0 

(96) 


we  see  that  as  x  the  rightward  propagating  wave  p(x,(jJ)  corresponding  to 
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the  probing  waves  (95)  takes  the  form 

p(x,00)  =  FCu)  exp-jcox,,  (97) 

so  that  the  Cholesky  recursions  provide  an  approximation  of  the  spectral 
factor  F(u)  of  w(u) . 

The  inverse  scattering  formulation  of  the  linear  least-squares  estimation 
problem  that  we  have  described  above  can  also  be  used  to  study  the  properties 
of  orthogonal  ladder  filters  of  the  type  described  in  Fig.  6b.  For  example, 
the  stability  and  lack  of  sensitivity  of  these  filters  to  roundoff  errors  are 
a  direct  consequence  of  the  losslessness  property  of  the  scattering  meditim 
[4] ,  [42] .  These  properties,  as  well  as  the  modularity  and  pipelinability  of 
ladder  filters  have  motivated  their  widespread  use  for  adaptive  equalization 
[43] ,  speech  processing  [12] ,  [44] ,  and  spectral  estimation  [14] ,  [45] . 

VI .  Inverse  Scattering  for  Asymmetric  Two-Component  Wave  Systems 

In  this  section,  the  inverse  scattering  problem  for  asymmetric  two- 
component  wave  equations  is  examined,  and  solved  by  using  two  coupled  sets 
of  Schur  recursions.  The  systems  which  are  described  by  asymmetric  two- 
component  wave  equations  are  not  necessarily  lossless,  and  we  can  therefore 
use  these  equations  to  describe  a  larger  class  of  physical  phenomena  than 
those  that  we  have  studied  in  the  previous  sections.  Our  results  will  be 
illustrated  by  considering  the  inverse  problem  for  a  nonuniform  transmission 
line  with  losses.  It  is  worth  noting  that  a  solution  of  the  inverse 
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scattering  problem  for  asymmetric  two-component  wave  equations  was  presented 
in  [10]  and  was  used  by  Jaulent  [46]  to  solve  the  inverse  problem  for  lossy 
transmission  lines.  However,  this  method  relied  on  the  solution  of  two 
coupled  Marchenko  equations,  whereas  the  solution  that  we  present  here  is 
differential,  and  uses  the  layer  stripping  principle  that  we  have  been 
advocating  throughout  this  paper. 

The  system  that  we  consider  is  described  by  the  asymmetric  two-component 
wave  equations 


which, 


d 

dx 


-jti)  -s  (x) 
-r(.x)  jO) 


P 

/s. 

q 


in  the  time-domain  correspond  to 


p^  +  p^  =  -s(x)  q(x,t) 
q^  -  q,  =  -r(x)  p{x,t) 


(98) 


(99a) 

(99b) 


It  is  assumed  that  r(x)  =  s (x)  =  0  for  x  <  0,  and  that  r,  s  8  L^[0,“) , 

so  that  r(x)  and  s  (x)  are  localized,  i.e.  they  go  to  zero  as  x  . 

Then,  the  scattering  matrix  5(03)  can  be  defined  as  in  Section  II  by 

relating  the  outgoing  and  incoming  waves  appearing  in  Fig.  2.  In  addition, 

T 

the  property  (9)  for  the  VJronskian  of  two  independent  solutions  a^(x,a3)  = 

(p.  (x,0))  ,  q.  (x,a3))  i=l,  2  of  (98)  remains  valid,  and  by  applying  it  to  the 
1  1 

waves  a^(x,03)  and  appearing  in  Figs.  2a  and  2b  respectively,  we 

obtain  the  reciprocity  relation 

(03)  =  f  (03) 

Li  K 


(100) 
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T 

However,  if  a  (x,a))  =  (p(!x,aj)_,  q(x,a)))  is  an  arbitrary  solution  of  (98), 
we  have 

^  (  |p|^  ■  1^1^)  =  2(.r;(:x)  -  s  (x) )  Re  (;^(x,u)q*  (x,a)) )  (101) 

so  that  the  scattering  mediiam  associated  to  (98)  is  not  lossless  unless 
r(x)  =  s (x)  which  corresponds  to  the  case  when  the  two-component  wave 
equations  are  symmetric.  This  implies  that  S(a))  is  not  a  unitary  matrix, 
and  consequently  we  cannot  recover  S  (co)  from  the  knowledge  of  the  left 
reflection  coefficient  only. 

Inverse  Scattering  Procedure 

The  inverse  scattering  method  that  we  develop  here  relies  on  the 
observation  that  if  time  is  reversed  (i.e.  t  is  changed  to  -t  in  (99)  ,  or 
0)  is  changed  to  -u  in  (98) ). ,  and  if  the  waves  P  and  q  are  interchanged, 
we  obtain  an  asymmetric  two-component  wave  system 

p^  +  p^  =  -r(x)  q^(x,t)  (102a) 

q^  -  q^  =  -s (x)  p^(x,t)  (102b) 

where  r(x)  replaces  s (x)  and  vice-versa.  The  scattering  matrix  associated 
to  this  system  is 


~0  1 

"o  l” 

A 

-1 , 

s  m  = 

_1  0. 

s  (-0)) 

1  0 

(S^(W)) 


(103) 
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where  to  obtain  (103)  we  have  used  the  reciprocity  relation  (100) .  The 
system  (102)  is  a  fake  system,  which  does  not  exist  really,  but  its  scattering 
matrix  is  entirely  specified  by  the  knowledge  of  S  (OJ)  . 

Then,  in  order  to  reconstruct  r(x)  and  s (x) ,  we  assume  that  the  true 
system  (98)  and  the  fake  system  (102)  are  probed  simultaneous ly  by  some 
waves  which  have  the  form 


p(x,t)  =  6 (t-x)  +  p(x,t)  u(t-x) 

q(x,t)  =  q(x,t)  u(t-x) 


(104) 


and 

p^(x,t)  =  5  (t-x)  +  p'^(x,t) u(t-x) 
A  ~A 

q  (x,t)  =  q  (x,t)u(t-x) 


(105) 


By  substituting  these  waves  in  (.98)  and  (102)  ,  we  obtain  the  system  of  coupled 
fast  Cholesky  recursions 


p  +  p  =  -s(x)  q(x,t) 

X  T- 


q^  -  q^  =  -r(x)  p(x,t) 


(106a) 


and 


with 


~A  -A  /  V  -A ,  ^ . 

p  +  p  =  -r(x)  q  (x,t) 

X  u. 

gA  ~A  ,  ,  ~A 

-  q.  =  -s(x)  p  (x,t) 

X  "C 

--  ~A 

r(x)  =  2q(x,x) ,  s (x)  =  2q  (x,x) 


(106b) 


(106c) 
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which  can  be  propagated  recursively  for  increasing  values  of  x,  starting 
from  X  =  0.  The  specification  of  the  initial  conditions  for  these  recursions 
is  very  important,  since  as  noted  above,  the  system  (102)  does  not  exist 
really  and  cannot  be  relied  upon  to  provide  some  experimental  waves  p^(0,t) 
and  q  (0,t) . 

The  initial  conditions  that  we  select  are 

p(0,t)  =  p^(0,t)  =  0  (107a) 

q(0,t)  =  R^(t)  ,  q^(0,t)  =  R^(t)  (107b) 

A 

where  R  (t)  and  R  (t)  denote  the  inverse  Fourier  transforms  of  the  left 

Li  J-i 

reflection  coefficients  R  (co)  and  R  (to)  -R  (to)  can  be  measured  directly, 

i-i  Ll  Xi 

and  from  (103) 

fij^((0)  =  (S~”(00))2^  (108) 

i.e.  :^(to)  is  the  (2,  1)  entry  of  the  inverse  of  S^(to)  .  Thus,  R^ (to)  can  be 
ll  L 

expressed  as  a  function  of  the  whole  scattering  matrix  S  (to)  ,  and  it  will  be 
specified  provided  that  we  can  measure  all  the  entries  of  S  (to)  .  This  implies 
that  we  must  have  access  to  both  ends  of  the  scattering  medium.  In  some  cases, 
such  as  for  the  inverse  problem  of  geophysics,  this  is  impossible;  but  for 
some  other  problems,  such  as  for  the  reconstruction  of  nonuniform  transmission 
lines ,  the  medirun  can  be  probed  from  both  sides ,  and  all  the  entries  of  S  (to) 
can  be  measured. 

Instead  of  expressing  our  reconstruction  procedure  in  terms  of  the 
coupled  fast  Cholesky  recursions  described  above,  we  can  use  a  set  of  coupled 
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Schur  recursions .  Let 


R  (.X ,  u) 


q(x,h3) 

p(x,a3) 


and 


(x,w) 


q  (x,(iO) 
pA(x,C0) 


(109) 


be  the  left  reflection  coefficients  for  the  true  and  fake  systems  over  the 
interval  [x,“)  ,  where  the  waves  p,  4',  p^,  in  the  definition  (109)  are 
assimied  to  have  the  form  (104)  -  (105)  .  Then,  R(x,u)  and  R  (x,(ja)  satisfy 
the  Riccati  equations 

R^  =  2j(JJR  +  s(x)  R  -  r(x)  (110a) 

“^A  •^A2 

R^  =  2ju)R  +  r(x)  R  -  s  (x)  (110b) 

with  intial  conditions 

R(0,w)  =  R_  (00)  ,  R^(0,a))  =  Rf(oo)  .  (Ill) 

L  Xj 

By  using  the  initial  value  theorem  for  the  reflection  coefficient  (109) ,  and 
taking  into  account  the  form  of  the  waves  (104)  -  (105) ,  we  get 

lim  2jcoR(x,co)  =  r(x)  (112a) 

lim  2jooR'^(x,oo)  =  s  (x)  (112b) 

(jVX>° 

which  can  be  combined  with  (110a)  and  (110b)  to  propagate  R(x,0J)  and  R  (x,oo) 
recursively,  and  to  reconstruct  r(x)  and  s (x)  for  all  x.  This  algorithm 
constitutes  the  generalization  of  the  Schur  algorithm  (21)  -  (22) . 

Reconstruction  of  Non-uniform  Transmission  Lines  with  Losses 

In  Section  III,  the  reconstruction  problem  for  a  nonuniform  lossless 


transmission  line  was  solved  using  the  Schur  algorithm.  We  now  consider 
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the  more  general  case  where  some  losses,  in  the  form  of  series  and  shunt 
resistances  per  unit  length  have  been  addded  to  the  transmission  line. 

This  reconstruction  problem  is  then  solved  as  an  asymmetric  two-component 
inverse  scattering  problem,  using  the  method  obtained  at  the  beginning  of 
this  section.  The  problem  is  set  up  as  in  [46]. 

An  infinitesimal  section  of  the  line  is  shown  in  Fig.  7.  R(x)  is  the 
non-uniform  series  resistance  per  unit  length,  representing  the  finite 
resistance  of  the  wires,  and  G(x)  is  the  shunt  conductance  per  unit  length, 
representing  leaikage  current  between  the  wires.  The  circuit  equations  are 


v(x,t)  =  (Li^  +  Ri)A  +  v(x+A,  t) 
i(x,t)  =  (Cv^  +  Gv)A  +  i(x+A,  t) 


(113) 


Dividing  by  A,  and  letting  A  -»■  0  yields  the  transmission  line  equations 


V  +  Li^  +  Ri  =  0 
X  t 


i  +  Cv  +  Gy  =0 

X  t 


(114) 


As  in  Section  III,  we  replace  the  position  x  by  the  travel  time  z (x) 
given  by  (39) ,  and  we  introduce  the  dimensionally  equivalent  variables 


V(z,t)  =  2  ^'^^v(z,t),  I(z,t)  =  Z^/^i(x,t) 


,1/2 


(115) 


where  Z(z)  =  (L(z)/C(z))  is  the  characteristic  impedance.  Then,  the 
equations  (114)  take  the  form 

R 

1  -  m  iz;  V 

(116) 


V  +1=-^  I  -  m(z)V 
Z  t  L 


I  +  =  m(z)  I  -  ^  V 

z  t  C 


where 


1  A. 

2  dz 


£n  Z ( z ) 


m(z) 


(117) 
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Making  the  change  of  variahles 


p(z,t)  =  j  (V+I),  q(z,t)  =  i  (V-I) 


(118) 


gives 


Pz  +  Pf  =  -a(z)p(z,t)  -  (m(z)  +  b(z))q{z,t) 


q_  -  q.  =  -  (m(z)  -  b(z))p(z,t)  +  a(z)q(z,t) 


(119a) 

(119b) 


which  is  almost  in  the  desired  form,  and  where 


,  ,  1  ,G  ^  R,  ,  ,  ,  1  ,G  R. 

a(z)  =  2  =  2 


(120) 


Considering  the  scaled  variables 


p^(z,t)  =  p(z,t)  exp 


rz 


a  (u)  du 


(121a) 


q^(z,t)  =  q(z,t)  exp  - 


a  (u)  du  , 


0 


(121b) 


and  taking  Fourier  transforms  yields  the  asymmetric  two-component  wave 
equations 


where 


1  /«3, 

p^  =  -jwp  (.z,w)  -  s(z)q  (z,a)) 


q^  =  -  r(z)p^(z,co)  +  joa  q^(z,a3) 


(122) 


f  z 


r(z)  =  (m-b)  exp  -2 


a(u)du  =  (^^(^n|) 


1  ,G  R, , 

2 


(|+  |)du 


(123a) 


rz 


s(z)  =  (m+b)  exp  2 


/  \  -j  ,1  d  ,  rt  hx  1  ,G  R.  . 
a(u)du  =  (4  ^  (^n  -)  +-(---))  exp 


. R  Gv  _ 

(-  + j)au 


(123b) 


-44- 


Thus,  if  we  are  given  the  scattering  matrix  S(w)  associated  to  the 
system  (122),  the  coupled  fast  Cholesky  or  Schur  recursions  (106)  and 
(110)- (112)  may  be  used  to  reconstruct  the  rather  bizarre  quantities  r(z) 
and  s(z).  Further,  these  two  quantities  are  the  most  information  about 
the  line  that  can  be  obtained  from  this  data.  Although  r(z)  and  s(z)  may 
seem  to  be  peculiar  quantities,  this  result  is  in  agreement  with  [46]. 

Note  that  in  the  event 

R(z)/L(z)  =  G(z)/C(z)  (124) 

we  may  recover  Z(z)  and  R(z)/L(z)  by  multiplying  and  dividing  r(z)  and  s(z), 

and  then  solving  two  differential  equations.  Thus,  in  this  case  it  is 

possible  to  recover  R(z),  L(z),  C(z) ,  and  G(z)  in  various  ratios  quite 

easily.  This  case  is  referred  to  as  the  Heaviside  condition  for  a  distortionless 

line  [32]  ,  since  if  (124)  holds  then  the  true  characteristic  impedance 

1/2 

( (R  +  jojL) /(G  +  jojC) )  which  relates  the  current  and  voltage  for  a  wave 
travelling  down  the  line,  is  real.  Thus,  the  current  and  voltage  for  such  a 
wave  are  in  phase,  just  as  in  the  lossless  line,  and  it  is  not  surprising 
that  ratios  of  various  line  parameters  can  be  recovered,  as  in  the  lossless 
case. 

VII .  Conclusion 

In  this  paper,  the  widespread  applicability  of  the  fast  Cholesky  and 
Schur  recursions  for  the  study  of  inverse  scattering  problems  has  been 
demonstrated.  These  algorithms  were  derived  by  using  a  layer  stripping 
principle  to  reconstruct  a  scattering  medium  described  by  symmetric  two- 
component  wave  equations,  for  the  case  when  the  medium  is  probed  by  impulsive 
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waves.  The  applicability  of  these  algorithms  to  the  reconstruction  of  a 
nonuniform  lossless  transmission  line ,  and  to  the  inverse  problem  for  a 
one-dimensional  layered  acoustic  medium  was  demonstrated.  In  addition, 
it  was  shown  that  the  linear  least-squares  estimation  problem  for  a 
stationary  process  could  be  posed  as  an  inverse  scattering  problem,  and 
solved  by  the  Schur  algorithm. 

Next,  an  asymmetric  two-component  inverse  scattering  problem  was 
considered  and  solved  by  using  a  coupled  set  of  fast  Cholesky  or  Schur 
recursions.  This  was  then  applied  to  the  inverse  problem  for  non-uniform 
transmission  lines  with  losses. 

There  are  several  topics  which  have  not  been  discussed  in  this 
paper  and  which  deserve  further  investigation.  One  of  them  is  the  study 
of  the  numerical  properties  of  the  Schur  algorithm  in  the  presence  of  noise 
or  modelling  uncertainties.  The  discrete-parameter  Schur  algorithm  was 
recently  shown  to  be  numerically  stable  by  Bultheel  [47] .  However,  a 
n\americally  stable  algorithm  can  perform  poorly  if  it  operates  on  ill 
conditioned  data,  which  could  happen  for  several  of  the  physical  problems 
that  we  have  examined  in  this  paper.  This  issue  deserves  therefore  to  be 
addressed.  An  additional  feature  of  the  layer  stripping  principle  that  we 
have  used  here  to  derive  the  fast  Cholesky  and  Schur  recursions  is  that  it 
is  quite  general,  and  it  is  applicable  to  more  general  physical  systems  than 
those  described  by  second  order  differential  equations.  For  example,  in 

[48],  ,  [49] ,  it  is  shown  that  this  principle  can  be  applied  to  the  reconstruction 
of  a  one-dimensional  elastic  medium  described  by  four  coupled  first-order 
differential  equations.  A  natural  extension  of  this  result  would  be  to  the 
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study  of  general  Hamiltonian  systems. 
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FIGURE  CAPTIONS 


Fig-  1: 

Fig.  2; 

Fig.  3: 
Fig.  4: 
Fig.  5: 
Fig .  6 : 


Elementary  scattering  sections  obtained  by  discretizing  the 
two- component  wave  equations. 

(a)  Scattering  for  an  impulsive  wave  incident  from  the  left, 
and  (b)  from  the  right. 

Infinitesimal  section  of  a  lossless  non-uniform  transmission  line. 

The  perceived  load  to  the  right  of  x. 

Inverse  Problem  for  a  layered  acoustic  medium. 

(a)  Aggregate  modeling  filter  for  y(-),  and  (b)  infinitesimal 
ladder  sections  associated  to  the  Krein- Levinson  recursions . 


Fig.  7:  Infinitesimal  section  of  a  lossy  non-'uniform  transmission  line. 
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